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Abstract 

The harmonic Frenkel-Kontorova model is used to ihustrate with an exactly solvable example 
a general technique of mapping a coherently strained epitaxial system with continuous atomic 
displacements onto a lattice gas model (LGM) with only discrete variables. The misfit strain of 
the original model is transformed into cluster interatomic interactions of the LGM. The clusters are 
contiguous atomic chains of all lengths but the interaction strength for long chains is exponentially 
small. This makes possible the application of efficient Monte Carlo techniques developed for discrete 
variables both in kinetic and equilibrium simulations. The formalism developed can be applied to 
ID as well as to 2D systems. As an illustrative example we consider the problem of self-organization 
of ID size calibrated clusters on the steps of the vicinal surfaces. 

PACS numbers: 05.65.+b, 68.66.La, 81.16.Dn 
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I. INTRODUCTION 



The phenomena of self-assembly and self-organization of coherent (i. e., dislocation-free) 
size calibrated nano- and atomic-scale structures observed during the heteroepitaxial growth 
in some systems§'i are considered to be promising tools for fabrication of microelectronic 
devicesi. 

A major factor influencing the phenomenon of self-assembly is the lattice size misfit 
(LSM) between the substrate and the growing overlayer commonly encountered in het- 
eroepitaxial systemsi. The LSM-induced strain is believed to be the driving force behind 
the above phenomena^'!. So an adequate account of the strain effects should lie at the 
basis of any theory of strained epitaxy. Because strained systems exhibit complicated ki- 
netics and morphologies, analytic approach is difficult, so a major technique in theoretical 
studies of strained epitaxy is the kinetic Monte Carlo simulation. The application of this 
technique, however, is severely hampered by the necessity to simulate the continuous atomic 
displacements. Because of this, atomistic models in such simulations are currently restricted 
to rather small systems consisting of only a few thousand atoms! while experimentally ob- 
served 3D quantum dots sometimes consist of several tens of thousand atoms eachi. 

Our research is based on the observation that as long as we are interested only in the 
coherent structures, there is a possibility to map the system onto a purely lattice model 
because in the absence of dislocations the relaxed atom can only either deviate from its 
symmetric position inside the same cell or to be displaced to another cell but there always 
exists a lattice site of a regular lattice to which this atom can be ascribed. So our first goal 
is to develop a formalism which would allow to map a coherent heteroepitaxial system with 
the continuous variables onto a lattice gas model with only discrete variables. 

One of important goals of the heteroepitaxial studies is the development of techniques of 
growth of ID quantum wires (QW) which can be used, e. g., for experimental investigation 
of the Luttinger model of interacting ID electronJ^. Furthermore, the QW may find appli- 
cation in microelectronics circuitryi, as well as in magnetic memory devicesi. The latter 
applications would require the ID structures of finite length. In the case of the memory 
devices it is also desirable that these atomic clusters (or chains) were of similar length and 
that they were arranged into periodic arrays in order to simplify the memory access. All 
these requirements can be satisfied by self-organized size calibrated structures similar to 
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quantum dots of the 2D epitaxyo. A phenomenological theory explaining the mechanism of 
formation of quantum dots was proposed in Ref.i. The theory is quite general and can be 
applied to objects in any number of dimensions. So to illustrate the techniques developed 
in the present paper we will study the conditions of formation of the ID size calibrated 
monatomic chains. 

II. THE MODEL 

A simple approach to theoretical description of the misfit is provided by the ID 
Frenkel-Kontorova model (FKM) which is frequently been used in qualitativei'0 and semi- 
quantitative studied! of strained epitaxy. To illustrate the generality of our approach we 
first present the formalism appropriate to the 2D epitaxy but in concrete examples we will 
restrict ourselves to the ID case. 

We consider an ensemble of a fixed number (A^) of atoms coherently deposited on a surface 
with a square lattice of deposition sites. The rectangular lattice geometry was chosen because 
it allows for the separation of x and y variables (see Ref.i and below). Let us first consider 
more general model with the energy of the epilayer of the form 



where Ri = {a^ix,a^iy) {a'^ the lattice constants in the two directions, ix and iy integers ), 
ni = 0, 1 is the occupation number of site i, Ui the atomic displacement, Vg the potential of 
interaction with the substrate, and V the interatomic potential. 

The analysis of the system considerably simplifies at low temperatures where it can be 
approximately reduced to a lattice gas model. This is achieved by exploiting the fact that 
the residence time of atoms at the lattice sites can be arbitrarily large due to the Arrhenius 
law obeyed by the probability of activated hopping over the potential barriers separating 
neighboring siteJii. The dynamics of the variables Ui, on the other hand, do not have any 
energy barriers. So at sufficiently low temperature these variables are capable of reaching 
their thermal equilibrium distribution during the time intervals between the atomic hops, 
i. e., with the atomic configuration remaining unchanged. Averaging over Ui will leave us 
with an effective non-equilibrium free energy function F^s of variables Ui only: 





{ni=l} 



n 




(1) 
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This purely lattice model can be further used in both equilibrium and kinetic studies. 

For the purposes of qualitative analysis it will suffice to average out the displacement 
variables in Eq. (|1|) in the harmonic approximation, i. e., by expanding Vs and Vp in the 
above equation up to the second order in the displacement variables UiS'El which is valid 
for small \u]\/a'^. This allows to perform the calculation exactly (see below). However, 
the integration can be extended to approximately account also for anharmonic terms. This 
extension, besides making the approach more accurate, can also account for some qualitative 
phenomenal. 

To facilitate comparison with other studies based on the Frenkel-Kontorova modeli'000 
we write the harmonic approximation as the second order power series expansion in Ui for 
Vs and in Ui — Ui_^ for the pair potential Vp-. 

+ \ '^pii'^i - «i+7 ~ A)^ ~ /7>i«i+7 
i7 

= Kiv + ^ E ^j^i^j 

ij 

+ 1 E k;Dlu7u] + Y: klMVln^n,^, (2) 

ij7 17 

where Vg and Vij are the values of the substrate and the pair potentials at the lattice sites 
for the atom in the symmetric position (the zero order approximation), 7 denotes x or 
y component; = —V!^/V!l^ (where V is the pair potential) is interpreted as the misfit 
parameter in the direction 7, kj and k'^ are the second derivatives of the corresponding 
potentials, is the dimensionless dynamical matrix defined by this equality, and V^-^i = 
F\ — F\^^. A major simplification is achieved under the nearest neighbor (NN) approximation 
for the relaxation because, as is seen from the last line of Eq. (0), in this case the variables 

and u\ separate, and the relaxations along the two directions are independent. 

With approximation (^ the statistical average amounts to the Gaussian integration to 
give 

^ ij ^7 

k p 

— Y' E G'ij V>i^i+7V]'njnj+^, (3) 

ij7 
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where G"' = and in Vg we gathered all terms which are proportional to the total particle 
number, such as the normalization of the determinant coming from the in-plane Gaussian 
integration and a similar term from the integration along z-direction. We do not discuss 
these contributions here because in the present study we restrict ourselves to systems with a 
fixed number of deposited atoms A^. In case of necessity these terms can be easily recovered. 
We also note the entropic contribution (the second term in the first line) which naturally 
appears in our formalism and which was shown to be crucial for proper description of the 
processes of depositionEl as well as for the correct prediction of the shape of the atomic 
clusterJli. 

Because of the gradient factors, only the ends of contiguous chains of atoms contribute 
into the last term of Eq. leaving two matrix elements of in the sum over ij inside 
every chain: the diagonal one which we denote as Cf''^ {I the length of the chain; we omit 
the superscript 7 to simplify notation) and the matrix element G'['~^'* connecting the two 
ends of the chain. Furthermore, in the NN approximation the matrix D is block-diagonal 
because the atoms belonging to different chains do not couple. Therefore, the determinant 
factorizes and the relaxation part of the free energy which consists of the terms in Eq. 
containing G takes the form 

i^relax = - E { ^ det + fc.f [G^ - Cf-^ ^1 } 
chains ^ 

= E i-TSi + Wi)^Y: Eu (4) 

chains chains 

where the summation is over the chains consisting of two or more atoms. Besides, in the 
first term on the right hand side (r. h. s.) we subtracted the single atom entropy term 
by assuming that it is accounted for in Vs which we discarded to simplify notation. It is 
easy to show that -Freiax can be expanded into an infinite sum of multiatom interactions of 
the form Vin-jiij^^^ . . . where 7 is the unit vector in the direction 7. The expansion 

coefficients are given by the expression 

Vi = Ei- 2Ei^i + Ei^2 (5) 

valid for all / > 2 with Eq = Ei = Q. 
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III. THE RELAXATION FREE ENERGY 



Omitting the superscript 7 to simplify notation, the matrix Di for chain of length / is 
defined according to Eq. (0) as 



kp^{Di)ijUiUj = ksu\ + kp{ui - U2Y + ksul + . . . + kp{ui-i - u{f + ksu] 
Hence, the square matrix kpDi have the following structure: 



(6) 



kpDi 



kri + kg 



-kn 



up I i\,g I up 




























/ijr) 2 Aj'Q I s 
























" kp kp ~\~ kg 



V 



(7) 



In Appendix |^ we have shown that the matrix elements of = 1/Di entering the expression 
for the relaxation free energy can be calculated with the use of recursion formulas for 
the tridiagonal matrices as [see Eqs. ( AlO ) and ( All )] 



6-;°^ = 1/(1 + a -di) 
Gfli = detCz+i = Cf^^bi, 

where a = kg/kp] di and bi are generated by the recursion relations 

di+i = l/{2 + a - di) and k+i = di^ik. 
The recursion relations (13) at large / drive di and hi to the fixed point di,,hi,: 



(8) 



(9) 



d^ 



1 + 



0. 



a 



Ia + -<1 



(10) 
(11) 



This means that at large / Cf'^ saturates to a constant value while = det Gi ~ d\, 
[see Eqs. (§)]. This in particular means that the entropic contribution into Freiax which is 
proportional to ln(detG/) is asymptotically linear in I (see Fig. |l]). 



Two contributions into -Freiax calculated with these formulas are shown in Fig. The 
entropic contribution is negative because in our canonical ensemble formalism we discarded 
the entropy corresponding to isolated atoms which is reflected in the denominator Gi 
in the expression for the free energy. This means, in particular, that when kp —>■ 
the atomic displacements within the cell became mutually independent (the on-site pair 
interactions do not depend on the atomic positions inside the cell), so Gi — * Gil, where / is 
the unit matrix, and the entropy term tends to zero [see Eq. @ and Fig. ||] . When kp is not 
equal to zero, the atomic displacements became restricted by its neighbors, so their entropy 
diminishes. For large chains the interior atoms are practically in translationally invariant 
environment, so for long chains the entropy loss is proportional to /. The above restrictions 
on the atomic displacements are the more stringent the larger kp. This is reflected in the 
slope of Si which is larger for smaller a = kg/kp (see Fig. |l]). 

Because the entropic contribution in Eq. (H) is multiplied by T, at low temperature 
it can be neglected. The remaining term Wi is negative and saturates to a finite value 
at large I (see Fig. |l]b). From Eq. (|^) it follows that the NN interatomic interaction 
Vnn = + where the first term (V^n) is some "chemical" interaction and the second 
term is the unrelaxed repulsion due to LSM. Thus, the relaxation energy in Eq. (^ is the 
difference between the relaxed elastic energy of the chain Ef^i{l) and the unrelaxed one 
Eq{1) = {I — l)kpf'^ with both energies being positive, as it should be for elastic energies. 
Thus, 

Wi = EUl)-E^,^{l) 

and is negative because relaxation lowers the energy. Now, if is negative but such that 
Vnn is not too large, then the reduced chain energy [Wi + Vnn(^ — 1)]/^ will have a minimum 
corresponding to an equilibrium chain length. In two dimensions these chains will unify in 
square platelets. Thus, our model contains the mechanism for the size calibration of atomic 
clusters which was phenomenologically described in Ref.i. In 2D this conclusion as well as 
the formula of Ref.i were verified with the use of the Monte Carlo simulationlll. 

Because the entropic contribution Si (see Fig. ||a) is practically linear in / according to 
Eq. d^) which has the form of the discrete second derivative it essentially contributes only 
into the pair interaction V2 (the multiatom contributions were found to be ~ 10%). Thus, 
the relaxation entropy formally amounts to effective NN interatomic repulsion which grows 
linearly with temperature. The entropic forces of this kind were earlier discovered in alloys 
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FIG. 1: Relaxation entropy in units oi ks (a) and relaxation energy in units of kpf"^ (b) for chains 
of length /. At both figures the first curve from the horizontal zero axis corresponds to q = 2, the 
second to a = 0.5, and the third one to a = 0.2 (note difference in scale). 

(see Appendix E of Ref .0 and references to earlier literature therein) . 

At low temperature the main contributions into Vi for I > 3 come from Wi (Fig.|l|b). 
Because the r. h. s. of Eq. (|^) has the form of the discrete second derivative and the curve 
Wi{l) is concave, all multiatom interactions are repulsive. On the other hand, in the case 
under consideration when atoms are supposed to assemble into clusters, there should exist 
attractive interactions which in our model are necessarily the pair ones [see Eq. (^]. 

IV. EXAMPLES FROM METALLIC HETEROEPITAXIAL SYSTEMS 

To illustrate the above formalism with realistic examples of strained epitaxy we consider 
two metallic heteroepitaxial systems — Ag/Pt and Pt/Co — which currently are being actively 
studied both experiment allyElS0 and theoreticallyilil. For simplicity we considered ID 
case and used the geometry of Ref.il where the growth on the steps of the closed packed 
(111) vicinal surface was studied. The position of a deposited atom was relaxed to its 
equilibrium value in order to find the value kg as the second derivative of the potential 
near equilibrium. The many-body "potentials" and corresponding parameters were taken 
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TABLE I: Parameters corresponding to the Pt/Co heteroepitaxial system 



e 


a 


kpf (eV) 


Vnn (eV) 


0% 


0.32 


0.049 


-0.218 


2% 


0.25 


0.100 


-0.179 


3% 


0.22 


0.136 


-0.152 



from Ref.0 for the Ag/Pt system and from Ref.il for the Pt/Co system. These potentials 
were devised specifically to application in the heteroepitaxy and for the (111) surfaces so 
we expect our results below are quite reliable. Other parameters listed in Table |I^ were 
calculated for the atomic pairs relaxed only in vertical position because in our approach we 
need the first and second derivatives [see discussion after Eq. (Q)] calculated in the center of 
of symmetry of the ID cell. In both calculations only the NN interactions were taken into 
account because the NNN ones were found to be negligeable. 

According to our calculations the Ag/Pt system has the following parameters: (the energy 
unit is eV, the length unit A): V^^ ^ -0.57, V^j^^ ^ -8.8 ■ IQ-^, k^p ^ 0.72, and a ^ 3.7. 
The large value of a means that the relaxation of the strain is very weak (see Fig. 0) so 
there is no size calibration with the above parameters. The crucial parameter a, however, 
can strongly vary in different systems. For example, according to our estimates based on 
the potentials of Ref.ii, in Pt/Co a ~ 0.32 is an order of magnitude smaller. Because of 
this the system is quite close to the self-assembly but the misfit strain is still too small. To 
enhance the misfit, in our model calculations we assumed that the Co underlayer is further 
compressed (e. g., by means of deposition on an appropriate substrate) by the factor 1 — e. 
(see Table |^ and Fig. |^). But for e below the critical value ec ~ 2% there is no self-assembly 
at T = because there is no minimum in Ei(T = 0)/l (see Ref.i). However, because of the 
T-dependent entropic contribution in the effective energy (^), slightly below Cc the system 
does exhibit the self-assembly at intermediate temperatures which, however, disappears as 
T ^ (see Fig. |^). Thus, our model predicts a new phenomenon which may be called the 
entropy driven transient self-assembly. For higher values of strain the Ei{T = 0)/l curve 
does have a minimum (see Fig. ^) in which case the qualitative analysis of Ref.i fully applies. 
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FIG. 2: Length dependence of the reduced energy at zero temperature of the Pt monatomic chains 
for different values of compression e of the Co substrate. 

V. CONCLUSION 

In this paper we considered a simple model of strained ID epitaxy and have shown that 
similarly to the 2D case the mechanism of self-assembly and size calibration proposed in 
Ref.i is operative also in this case. We considered only two explicit examples for which 
there exist in literature reliable interatomic potentials and in one system already found the 
size calibration for stressed substrate. In our opinion, this shows that the above phenomena 
should be as common in ID as they are in 2D heteroepitaxy. Further argument in favor of this 
conclusion is that for the size calibration to be plausible the crucial parameter a = kg/kp 
should be as small as possible. This favors small values of kg- However, the geometry 
considered in the present paper is not quite favorable because of the high coordination 
of the atoms deposited at the steps with 5 NN atoms of the substrate (see which 
enhances kg. It may be hoped that with lower coordination of the deposited atoms (as in the 
case of ID structures of ReLQ) the conditions for the size calibration will be more favorable. 
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APPENDIX A 



According to Eq. 



D, 



( \ 

1 + a-l 0... 

-1 2 + a -1 ... 

-1 2 + a -1 ... 



... -1 2 + a -1 
... -1 2 + a -1 
... -1 1 + a 

\ 



V 



(Al) 



where a = kg/kp. Because matrices Di are tri-diagonal, their determinants satisfy recurrence 
relations which can be used to calculate all quantities entering Eqs. (^) and (^. The diagonal 
element of the matrix Gi = 1/ Di is 



GS°) = n_i/detA 



(A2) 



where r^^i is the determinant of the matrix obtained from Di by crossing out its first row 
and the first column: 



2 + a -1 ... 
-1 2 + a -1 ... 
-1 2 + a -10 



... -1 2 + a -1 
... -1 1 + a 



(A3) 
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Expanding det Di with respect to the first row we get 



det A = (1 + a)r;_i -n_2. (A4) 
Comparing this with Eq. ( [A^ ) we get 

= (A5) 

where di^i = ri^2/ri~i. Now, expanding determinant Eq. ([A3D with respect to the elements 
of the first row we get the following three term recurrence relation 

n_i = (2 + a)ri_2 - n_3. (A6) 



But because Eq. (K2) includes only the ratio ri^2/^i-i by dividing Eq. (|A6|) by ri^2 we can 



transform it to a simpler form 

= 2 + a - di^2- (A7) 

<li-i 

The off-diagonal matrix element Cf of the inverse matrix Gi is equal to the ratio of the 
determinant of the matrix obtained from Di by crossing out its first row and the last column 
multiplied by (—1)'+^ and divided by det/)/. As is easy to see from Eq. (|A1|) the matrix 
thus obtained is a triangular matrix with the diagonal elements being all equal to -1. Hence 
its determinant is equal to (—1)'"^ and 

Gf-^^ = 1/ det Di = det Gi. (A8) 

From Eq. (^5) we get 



G?-'^ = = = ^^^^ = k.,Gl'\ (A9) 

where = l/r^-i = di^ibi-2 (see the definition of di^i above). Finally, making the shift 
/ — 1 ^ / + 1 the above formulas can be summarized as 

= i/n + a-di) 

(AlO) 

G^li = det Gi+i = G^^^bi, 
where di and bi are generated by the recursion relations 

di+i = 1/(2 + a - di) and bi+i = di+ik (All) 
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initialized by do = ^'o = 1- 
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